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We introduce a generalized d-dimensional Fermi-Pasta-Ulam (FPU) model in presence of long- 
range interactions, and perform a first-principle study of its chaos for d = 1, 2, 3 through large-scale 
numerical simulations. The nonlinear interaction is assumed to decay algebraically as d“.“ (a > 0), 

{dij} being the distances between N oscillator sites. Starting from random initial conditions we 
compute the maximal Lyapunov exponent Xmax as a function of N. Our >> 1 results strongly 
indicate that Xmax remains constant and positive for a/d > 1 (implying strong chaos, mixing and 
ergodicity), and that it vanishes like N~'^ for 0 < a/d < 1 (thus approaching weak chaos and 
opening the possibility of breakdown of ergodicity). The suitably rescaled exponent k exhibits 
universal scaling, namely that (d -|- 2)k depends only on a/d and, when a/d increases from zero 
to unity, it monotonically decreases from unity to zero, remaining so for all a/d > 1. The value 
a/d = 1 can therefore be seen as a critical point separating the ergodic regime from the anomalous 
one, K playing a role analogous to that of an order parameter. This scaling law is consistent with 
Boltzmann-Gibbs statistics for a/d > 1, and possibly with g-statistics for 0 < a/d < 1. 

PACS numbers: 05.70.-a, 05.45.Pq, 05.45.-a, 89.75.Da 


I. INTRODUCTION 

Many-body systems with long-range-interacting forces 
are very important in nature, the primary example being 
gravitation. Long-ranged systems deviate significantly 
from the conventional ‘well behaved’ systems in many 
respects. Various features like ergodicity breakdown, 
ensemble inequivalence, non-mixing nonlinear dynam¬ 
ics, partial (possibly hierarchical) occupancy of phase 
space, thermodynamical nonextensivity for the total en¬ 
ergy, longstanding metastable states, phase transitions 
even in one dimension, and other anomalies, can be ob¬ 
served in systems with long-range interactions. Consis¬ 
tently, some of the usual premises of Boltzmann-Gibbs 
(BG) statistical mechanics are challenged and an alter¬ 
native thermostatistical description of these systems be¬ 
comes necessary in many instances. For some decades 
now, g-statistics [ 3,0 has been a useful formalism to 
study such systems, and has led to satisfactory exper¬ 
imental validations for a wide variety of complex sys¬ 
tems (see for instance M)- The deep understanding 
of the microscopical nonlinear dynamics of such systems 
naturally constitutes a must in order to theoretically le¬ 
gitimize the efficiency of the (/-generalization of the BG 
theory. For classical systems such as many-body Hamil¬ 
tonian ones and low-dimensional maps, a crucial aspect 
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concerns the sensitivity to the initial conditions, which is 
characterized by the spectrum of Lyapunov exponents. If 
the maximal Lyapunov exponent Xmax is positive, mixing 
and ergodicity are essentially warranted, and we conse¬ 
quently expect the BG entropy and statistical mechanics 
to be applicable. If instead Xmax vanishes, the sensitiv¬ 
ity to the initial conditions is subexponential, typically a 
power-law with time, and we might expect nonadditive 
entropies such as Sq and its associated statistical me¬ 
chanics to emerge, as has been observed numerically as 
well as experimentally in many systems (see, for instance, 

iiiil). 


II. MODEL AND THE NUMERICAL SCHEME 


In the present paper we extend to d-dimensions {d = 
1,2,3) and numerically study from first principles (i.e., 
using only Newton’s law F = ma) the celebrated Fermi- 
Pasta-Ulam (FPU) model with periodic boundary condi¬ 
tions; nonlinear long-range interactions between all the 
N = L‘^ oscillators are allowed as well. The Hamiltonian 
is the following one: 
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( 1 ) 

where r/ and pi are the displacement and momentum 
of the i-th particle with mass rm = m; a > 0, b > 0, 
and a > 0. Here dij is the shortest Euclidean distance 
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between the i-th and j-th lattice sites (1 < i,j < N)-, this 
distance depends on the geometry of the lattice (ring, 
periodic square or cubic lattices). Thus for d = 1, dij = 
1,2,3,...; for d = 2, dij = 1,-^2, 2,..., and, for d = 3, 
dij = 1, V2, 73, 2,... If a/d > 1 (0 < a/d < 1) we have 
short-range (long-range) interactions in the sense that 
the potential energy per particle converges (diverges) in 
the thermodynamic limit N ^ oo; in particular, the a —>• 
c» limit corresponds to only first-neighbor interactions, 
and the a = 0 value corresponds to typical mean field 
approaches, when the coupling constant is assumed to be 
independent from distance. The instance (d, a) = (1,cxd) 
recovers the original /3-FPU Hamiltonian, that has been 
profusely studied in the literature; the d = 1 model and 
generic a has been addressed i n 1241 . 

Although not necessary (see [1^), we have followed the 
current use and have made the Hamiltonian extensive for 
all values of a/d by adopting the scaling factor N in the 
quartic coupling, where 

( 2 ) 

i=i “b 

hence, iV depends on a,N,d, and the geometry of the 
lattice. Note that for a = 0 we have N = N, which 
recovers the rescaling usually introduced in mean field 
approaches. In the thermodynamic limit N —>■ c», N 
remains constant for a/d > 1, whereas N ~ ^-a/d 

0 < a/d < 1 {N ^ InA^ for a/d = I); see details in [l^ 
and references therein. 

Let us mention that the analytical thermostatistical 
approach of the present model is in some sense even 
harder than that of coupled XY or Heisenberg rotators 
already addressed in [l^ 1251 - 1^ . Indeed, the standard 
BG approach of these models is analytically tractable, 
whereas not even that appears to be possible for the 
original FPU, not to say anything for the present gener¬ 
alization. Therefore, for this kind of many-body Hamil¬ 
tonians, the numerical approach appears to be the only 
tractable one. 

To numerically solve the equations of motion (New¬ 
ton’s law) we have employed the symplectic second or¬ 
der accurate velocity Verlet algorithm. To accelerate the 
computationally expensive part of the force calculation 
routine we have exploited the convolution theorem and 
used a Fast Fourier transform algorithm. This yields a 
considerable reduction in the number of operations for 
force calculation from 0{N‘^) to 0{NInN), thus facil¬ 
itating computation for larger system sizes and longer 
times. 

We choose the time step At (which is typically ~ 10“^ 
for most of our results) such that the standard devi¬ 
ation of the energy density over the entire simulation 
time (i.e., the number of iterations required by the max¬ 
imal Lyapunov exponent to saturate, which is typically 
^ 10® — 10® iterations, depending on system parameters) 
is of the order of 10“"^ or smaller (for the range of N 
considered here, 10 < N < 10®). 


Starting from a random initial displacements 7 drawn 
from a uniform distribution centered around zero, and 
momenta pi from a Gaussian distribution with zero mean 
and unit variance, we evolve the system and compute the 
maximal Lyapunov exponent Xmax defined as follows: 


Xmax = lim lim i In , 

t^ooS{0)^0t (5(0) 


( 3 ) 


where S{t) = is the metric distance 

between the fiducial orbit and the reference orbit having 
initial displacement <5(0). We numerically compute this 
quantity by using the algorithm by Benettin et al [^ . 
For typical values of the exponent a, we compute Xmax 
as a function of the system size N for d = 1,2, and 3. 


III. SIMULATION RESULTS 


Let us now present the results of our numerical anal¬ 
ysis by setting to = 1 (no loss of generality), and fixing 
the energy density u=U/N = 9.0 and b = 10.0 for all d, 
unless stated otherwise, where U is the total energy as¬ 
sociated with TL. Additionally, we have set the harmonic 
term to zero, i.e. a = 0, for reasons that will be elabo¬ 
rated later. In fact such a model, with only the quartic 
anharmonic nearest neighbor interactions, has been stud¬ 
ied previously in the context of heat conduction . 
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FIG. 1. (Color online) Log-log plot of the dependence, for 
a ring {d = 1), of the maximal Lyapunov exponent Xmax 
on the number N — L of oscillators for (a, b, u) = (0,10, 9) 
and typical values of the exponent a. Each individual curve 
has been multiplied by the number indicated next to it for 
visualization clarity. 


In Figs, m [2] and|3]we present, for d =1, 2 and 3 re¬ 
spectively, the maximal Lyapunov exponent Xmax as a 
function of the system size for typical values of the ex¬ 
ponent a. We find that, for a > d, Xmax saturates to a 
positive value with increasing N, which strongly suggests 
that it will remain so for A —oo, thus leading to ergod- 
icity, which in turn legitimizes the BG thermostatistical 
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theory. In contrast, for 0 < a < d, \max algebraically 
decays with N = L‘^ as 

Xmax ~ N-'^ (4) 

where k > 0 and depends on (a,d). Assnming 
that it remains so for increasingly large iV, we expect 
limjv^oo Xmax = 0, which implies that the entire Lya¬ 
punov spectrum vanishes. This characterizes weak chaos 
for 0 < a/d < 1, i.e., subexponential sensitivity to 
the initial conditions, which opens the door for break¬ 
down of mixing, or of ergodicity, or some other nonlin¬ 
ear dynamical anomaly. Within this scenario, the vio¬ 
lation of Boltzmann-Gibbs statistical mechanics in the 
TV —>■ oo limit becomes strongly plausible (see, for exam¬ 
ple, l2l,[2l). 
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FIG. 2. (Color online) The same as in Fig. 1 for a periodic 
square lattice {d = 2) with N = Lp’ oscillators. 



FIG. 3. (Color online) The same as in Fig. 1 for a periodic 
cubic lattice (d = 3) with N = oscillators. 

From the results illustrated in Figs. (TJ [2] and [3] we 
compute the exponent K(a, d) for d = 1, 2 and 3, as shown 
in Fig. m including its inset. We find that K(a, d) > 0 for 


0 < a < d, and, within numerical accuracy, vanishes for 
a > d. Also note that K(0,d) decreases for increasing d. 
Remarkably enough, all three curves in the inset of Fig. 
0 ] can be made to collapse onto a single curve through 
the scalings a —t a/d and K(a,d) — >■ (d -I- 2) /t(a, d). 
This is shown in the main figure of Fig. 0 ] In other 
words, (d + 2) K(a,d) = f{a/d) where f{x) appears to be 
a universal function. 

A similar scaling was also verified for the classical 
model of long-ranged coupled rotators [l^, . Some rel¬ 

evant differences exist however between the two models 
and their sensitivities to initial conditions. The long- 
range-interacting planar rotator model exhibits, for a 
critical energy density Uc [H, a second order 

phase transition from a clustered phase (ferromagnetic) 
to a homogeneous one (paramagnetic). Such critical phe¬ 
nomenon does not exist in either the short-ranged or 
the long-ranged FPU model. For the XY ferromagnetic 
model the exponent k, for a = 0 is found to be indepen¬ 
dent from d (quite obvious since the a = 0 model has no 
dimension) and given by /v(0,d) = 1/3 [1^ [3l| (see also 
[^). In contrast, our long-range model yields a value 
K(0,d) which depends on d. Indeed, for d = 1, 2 and 3, 
we respectively obtain k( 0, d) ~ 1/3, 1/4 and 1/5. 

This difference in k( 0, d) is related to the fact that, for 
the XY model, the number of degrees of freedom (num¬ 
ber of independent variables needed to specify the state 
of the system in phase space) for N coupled rotators in 
d dimensions is 2N (Vd), whereas, for our model, there 
are 2Nd degrees of freedom, hence the dimension of the 
full phase space grows linearly with d. Thus there are 
more possible phase space dimensions for our coupled os¬ 
cillator system to escape even if gets somewhat trapped 
in some non-chaotic region of the phase space. Conse¬ 
quently, the system gets closer to ergodicity (equivalently, 
K gets closer to zero) for increasing d. It is even not ex¬ 
cluded that, because of some generic reason of this kind, 
k( 0, d) (Vd) for the long-ranged XY model and k( 0, 1) for 
the system studied here, we obtain (in absence of the 
integrable term, i.e., with a = 0) the same value 1/3. 

In this context we should mention another recent study 
[ 3 ^ of the Hamiltonian mean field (HMF) model which 
is the a = 0 particular case of the long-ranged XY model 
discussed above. Using numerical and analytical argu¬ 
ments it was suggested that the nature of chaos is quite 
different for this model (which has a phase transition at 
Uc = 3/4) in the homogeneous phase {u > Uc) where 
Xmax the ordered phase {u < Uc) where Xmax 

remains positive and finite, and at criticality (it —>■ Uc) 
where Xmax ~ in the infinite size limit. How¬ 

ever in another earlier work [^ . using scaling argu¬ 
ments and numerical simulations, it was observed that 
Xmax below the critical point {u = 0.69) in the 

(non equilibrium) quasi-stationary regime of the HMF 
system. 

Another class of models might also have a similar be¬ 
havior. If we consider the d-dimensional long-range¬ 
interacting n-vector ferromagnet, we expect an expo- 
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nent K.{a,n,d). We know that for n = 2 (XY symme¬ 
try) k( 0, 2,1) = 1/3, for n = 3 (classical Heisenberg 
model symmetry) k( 0, 3,1) = 0.225 ± 0.030 [s^, and 
for n —^ oo (spherical model symmetry) most plausibly 
K(0,oo,d) = 0 (yd). These expressions can be simply 
unified through K{0,n,d) = l/(n -|- 1) (Vd). 

Strikingly enough, the present Fig. [Hand Fig. 2 of 
for the d-dimensional XY model are numerically indistin¬ 
guishable within error bars. This suggests the following 
heuristic expression: 


k(q;, d) 
k( 0, d) 


/(a/d) 


1 — (a/d)^ 

1 -I- (a/d)2/6 ’ 


( 5 ) 


where this specific analytic expression for f{x) has been 
first suggested in [2^. This or a similar universal be¬ 
havior is expected to hold for d-dimensional long-range¬ 
interacting many-body models such as the present one, 
the XY ferromagnetic one, and others such as, for in¬ 
stance, the n-vector ferromagnetic one (Vn). 
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FIG. 4. (Color online) The inset shows the exponent ii{a, d) as 
a function of a for d = 1, 2, 3. Note that k > 0 for 0 < a < d 
and K = 0 for a > d. The main figure exhibits the universal 
law obtained by appropriately rescaling the abscissas and or¬ 
dinates as indicated on the axes, i.e., (d-|-2)K(a, d) = /(a/d). 
The thick continuous curve is the heuristic scaling function 
f{x) = (1 — x^)/{l + a;^/6) [2^, which, within the present 
precision, is a remarkably close fit to the collapsed data. The 
present collapse obviously implies k(0, d) = l/(d-|- 2), hence 
limd_>oo «:(0, d) = 0, thus recovering ergodicity, as intuitively 
expected. 

All the numerical results presented until now are with 
a fixed set of parameters (a, 6, it,) = (0,10,9) and a fixed 
time step At. Before concluding, let us briefly mention 
some results concerning the influence of these parameters 
on Xmax{N) and K{a,d). In Fig. [5ji we plot Xmax{N) 
for d = 1 for three different sets of {b, it) keeping all 
other parameters unchanged. We find that increasing 
b has the same effect as increasing u ~ the maximum 
Lyapunov exponent Xmax increases with both of them 
but the slope of the curve k, remains practically unaltered. 
For a = 0, it is straightforward to show that the average 




of Hamiltonian Eq. o remains invariant with respect to 
b and u (all other parameters remaining the same) under 
the transformations 


x'= (b/u)^^'^ X, t'= {bu)^^^ t. (6) 

The second transformation in Eq. ([5]) implies that the 
maximum Lyapunov exponent Xmax ("^ ^~^) satisfies the 
following scaling relation: 


>^'max = >^max ■ (7) 

Using the data in the main figure, we show in the inset 
of Fig. [S^ the variation of X'max = {bu)~^('^ Xmax with 
N. As predicted by the scaling analysis, we get an excel¬ 
lent data collapse of the three curves. This is precisely 
as desired, keeping in mind the universal behavior ubiq¬ 
uitously found in statistical mechanics, in the sense that 
scaling indices, such as n here, are generically expected to 
be independent of the microscopic details of the model. 

For nonzero values of a, the simple scaling Eq. © dis¬ 
appears, and Xmax{N) shows a saturation to a positive 
value that vanishes for a = 0 when N is large, b being 
a finite positive number. This is shown in Fig. [5]3 for 
two values of a with the same value of b. The satura¬ 
tion of Xmax for a > 0 needs careful study to be under¬ 
stood properly. In Fig. [5]: we have shown (for d = 2) 
that increasing At can also lead to a deviation from the 
Xmax ~ A”” behavior; this deviation is quite expected, 
and one should choose the time step judiciously. Note 
that the saturation behavior in Fig. [SJi is not due to 
finiteness of the time step. 


IV. SUMMARY AND DISCUSSIONS 

Summarizing, we have introduced a d-dimensional 
generalization of the celebrated Fermi-Pasta-Ulam 
model which allows for long-range nonlinear interac¬ 
tion between the oscillators, whose coupling constant 
decays as distance~°‘. We have then focused on the 
sensitivity to initial conditions, more precisely on the 
first-principle (based on Newton’s law) calculation of 
the maximal Lyapunov exponent Xmax as a function of 
the number N of oscillators using large-scale numerical 
simulations. Without the quadratic nearest neighbor 
interaction (i.e., a = 0), Amax(A) appears to asymptot¬ 
ically behave as N~'^ (with k > 0) for 0 < a/d < I, 
and approach a positive constant (i.e., k = 0 ) for 
a/d > 1 in the A —> oo thermodynamic limit. Our 
results provide strong indication that k only depends 
on {a,d), and does so in a universal manner, namely 
(2 -I- d)K{a,d) = f{a/d) for 0 < a/d < I, and k = 0 
for a/d > 1. This universal suppression of strong chaos 
is well approximated by a model-independent heuristic 
function f{x) ~ (I — x^)/(I -b x‘^/6), previously found 
[a M for the d-dimensional XY model of coupled 
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FIG. 5. (Color online) Parameter dependencies of Xmax{N) with N = L‘^-. (a) for different b’s and u’s with {a, a, At) = 
(0,0,0.002) - inset shows data collapse obtained by rescaling the y-axis of the main figure as Xmax', (b) for different 

a’s with (q, 6 , u, At) = (0,10, 9, 0.002); (c) for different At’s with (a, a, b, u) = (0,0,10, 9). 


rotators. Thus, in the thermodynamic limit, these 
systems (and plausibly others as well) have a sort of 
critical point at a/d = 1, which separates the ergodic 
a/d > 1 region (where the Boltzmann-Gibbs statistical 
mechanics is valid, and the stationary state distribution 
of velocities is the standard Maxwellian one), from the 
weakly chaotic 0 < a/d < 1 region with anomalous 
nonlinear dynamical behavior (where g-statistics might 
be expected to be valid, and the one-body distribution 
of velocities appears to be of the g-Gaussian form, 
consistently with preliminary results available in the 
literature ). The present universality scaling for 

K{a/d) enables the conjecture that the indices q of the 
distributions of velocities and of energies might exist 


and only depend on the ratio a/d. Naturally, all these 
observations need further and wider checking, which 
would be welcome. Work along this line is in progress. 

Acknowledgments: We are thankful to H. 

Ghristodoulidi and A. Ponno for sharing with us their 
observations concerning the role of the coupling constant 
a in the determination of k. We have also benefitted 
from fruitful discussions with L.J.L. Girto, G. Sicuro, 
P. Rapcan, T. Bountis, and E.M.F. Curado. We grate¬ 
fully acknowledge partial financial support from GNPq 
and Faperj (Brazilian agencies) and the John Templeton 
Foundat ion-US A. 


[1] C. Tsallis, J. Stat. Phys. 52, 479 (1988); C. Tsallis, R.S, 
Mendes and A.R. Plastino, Physica A 261, 534 (1998). 

[2] M. Gell-Mann and C. Tsallis, eds., Nonextensive Entropy 
- Interdisciplinary Applications (Oxford University Press, 
New York, 2004); C. Tsallis, Introduction to Nonex¬ 
tensive Statistical Mechanics - Approaching a Complex 
World, (Springer, New York, 2009). 

[3] P. Douglas, S. Bergamini, and F. Renzoni, Phys. Rev. 
Lett. 96, 110601 (2006); E. Lutz and F. Renzoni, Nature 
Physics 9, 615 (2013). 

[4] B. Liu and J. Goree, Phys. Rev. Lett. 100, 055003 (2008). 

[5] R.M. Pickup, R. Cywinski, C. Pappas, B. Farago, and P 
Fouquet, Phys. Rev. Lett 102, 097202 (2009). 

[6] R.G. DeVoe, Phys. Rev. Lett. 102, 063001 (2009). 

[7] L.F. Burlaga, A.F. Vinas, N.F. Ness, and M.H. Acuna, 
Astrophys J. 644, L83 (2006). 

[8] G. Livadiotis and D.J. McGomas, Astrophys. J. 741, 88 

( 2011 ). 

[9] A. Upadhyaya, J.-P. Rieu, J.A. Glazier and Y. Sawada, 
Physica A 293, 549 (2001). 

[10] J.S. Andrade, Jr, G.F.T da Silva, A.A. Moreira, F.D. 
Nobre and E.M.F. Gurado, Phys. Rev. Lett 105, 260601 
( 2010 ). 

[11] CMS Gollaboration, Phys. Rev. Lett. 105, 022002 
(2010), and J. High Energy Phys. 8, 86 (2011). 


[12] ALIGE Collaboration, Phys. Lett. B 693, 53 (2010), and 
Phys. Rev. D 86, 112007 (2012). 

[13] ATLAS Collaboration, New J. Phys. 13, 053033 (2011). 

[14] PHENIX Collaboration, Phys. Rev. D 83, 052004 (2011), 
and Phys. Rev. C 84, 044902 (2011). 

[15] C.Y. Wong and G. Wilk, Phys. Rev. D 87, 114007 (2013). 

[16] L. Marques, E. Andrade-II and A. Deppman, Phys. Rev. 
D 87, 114022 (2013). 

[17] G. Combe, V. Richefeu, G. Viggiani, S. A. Hall, A. Ten- 
gattini and A. P. F. Atman, AIP Conf. Proc 1542, 453 
(2013); G. Gombe, V. Richefeu, M. Stasiak and A. P. F. 
Atman, Phys. Rev. Lett. 115, 238301 (2015). 

[18] M. L. Lyra and C. Tsallis, Phys. Rev. Lett. 80, 53 (1998). 

[19] G. Anteneodo and G. Tsallis, Phys. Rev. Lett. 80, 5313 
(1998). 

[20] F. Baldovin and A. Robledo, Phys. Rev. E 69, 045202(R) 
(2004). 

[21] G. Gasati, C. Tsallis and F. Baldovin, Europhys. Lett. 
72, 355 (2005). 

[22] E. P. Borges, C. Tsallis, G.F.J. Ananos and P.M.G. 
Oliveira, Phys. Rev. Lett. 89, 254103 (2002). 

[23] U. Tirnakli and E. P. Borges, Scientific Reports 6, 23644 
(2016). 

[24] H. Ghristodoulidi, C. Tsallis and T. Bountis, EPL 108, 
40006 (2014). 

[25] M. Antoni and S. Ruffo, Phys. Rev. E 52, 2361 (1995). 










6 


[26] A. Campa, A. Giansanti, D. Moroni, C. Tsallis, Phys. 
Lett. A 286, 251(2001). 

[27] L.J.L. Cirto, V.R.V. Assis and C. Tsallis, Physica A 393, 
286-296 (2014). 

[28] L. J. L. Cirto, L. S. Lima and F. D. Nobre, J. Stat. Mech. 
P04012 (2015). 

[29] G. Benettin, L. Galgani, J.-M. Strelcyn, Phys. Rev. A 
14, 2338 (1976). 

[30] S. Lepri, R. Livi, and A. Politi, Chaos 15, 015118 (2005). 


[31] M.C. Firpo, Phys. Rev. E 57, 6599 (1998); M.C. Firpo 
and S. Ruffo, J. Phys. A 34, L511 (2001); C. Anteneodo 
and R.O. Vallejos, Phys. Rev. E 65, 016210 (2001). 

[32] R. Bachelard and M. Kastner, Phys. Rev. Lett. 110, 
170603 (2013). 

[33] F. Ginelli, K. A. Takeuchi, H. Chate, A. Politi, and A. 
Torcini, Phys. Rev. E 84, 066211 (2011). 

[34] V. Latora, A. Rapisarda, C. Tsallis, Physica A 305, 129 

( 2002 ). 

[35] F.D. Nobre and C. Tsallis, Phys. Rev. E 68, 036115 
(2003). 


